# 学姐给的数据读取文件 读取、计算ETo数据的脚本

import pyeto
import os
import math

month_of_days31 = [1, 3, 5, 7, 8, 10, 12]
month_of_days30 = [4, 6, 9, 11]
feb_month = 2


def is_leap_year(year):
    """
    判断当前年份是不是闰年，年份公元后，且不是过大年份
    :param year: 年份
    :return: True 闰年， False 平年
    """
    if year % 4 == 0 and year % 100 != 0 or year % 400 == 0:
        return True
    return False


def get_day_of_year(year, month, day):
    """
    获取一个日期在这一年中的第几天
    :param year: 年份
    :param month: 月份
    :param day: 日期
    :return: 在这一年中的第几天
    """
    # 参数校验

    if month == 1:
        return day

    if month == 2:
        return day + 31

    days_of_31_num = 0
    days_of_30_num = 0
    # 31天月份数
    for days_of_31 in month_of_days31:
        if days_of_31 < month:
            days_of_31_num += 1
        else:
            break

    # 30天月份数
    for days_of_30 in month_of_days30:
        if days_of_30 < month:
            days_of_30_num += 1
        else:
            break

    return days_of_31_num * 31 + days_of_30_num * 30 + (29 if is_leap_year(year) else 28) + day


def pm(tmin, tmax, ws_10, tavg, latitude, altitude, year, month, day, rh_mean, sunshine_hours):
    """
    ETo 计算 / 日数据  estimating reference evapotransporation (ETo) for a grass reference crop
    :param tmin:  Temperature mini/day
    :param tmax: Temperature max/day
    :param ws_10: Measured wind speed [m s-1]   10m风速
    :param tavg: Temperature averate 平均气温
    :param latitude:
    :param altitude:
    :param year:
    :param month:
    :param day:
    :param rh_mean: Mean relative humidity [%] (average of RH min and RH max).
    :param sunshine_hours:  Sunshine duration [hours]
    :return: evapotranspiration (ETo) [mm day-1]  ；   Wind speed at 2 m above the surface [m s-1]
    """
    # tmin = -5.1
    # tmax = 1.3
    # tavg = -1.1  # 平均气温
    # ws_10 = 1.8  # 10m风速
    latitude = pyeto.deg2rad(latitude)  # 纬度
    # altitude = 31.3  # 海拔高度
    day_of_year = get_day_of_year(year, month, day)
    # rh_mean = 48  # 平均相对湿度
    atmos_pres = 101.85  # 大气压
    psy = 0.0677  # 湿度常数
    shf = 0.0
    coastal = False  # True or False 是否为沿海地区 不需要
    sol_dec = pyeto.sol_dec(day_of_year)  # 计算一年中的太阳赤纬
    sha = pyeto.sunset_hour_angle(latitude, sol_dec)
    ird = pyeto.inv_rel_dist_earth_sun(day_of_year)
    et_rad = pyeto.et_rad(latitude, sol_dec, sha, ird)
    cs_rad = pyeto.cs_rad(altitude, et_rad)
    daylight_hours = pyeto.daylight_hours(sha)

    sol_rad = pyeto.sol_rad_from_sun_hours(daylight_hours, sunshine_hours, et_rad)  # rs
    svp_tmin = pyeto.svp_from_t(tmin)  # 每日最低温度下的饱和蒸气压
    svp_tmax = pyeto.svp_from_t(tmax)  # 每日最高温度下的饱和蒸气压
    avp = pyeto.avp_from_rhmean(svp_tmin, svp_tmax, rh_mean)
    ni_sw_rad = pyeto.net_in_sol_rad(sol_rad, albedo=0.23)  # rns
    no_lw_rad = pyeto.net_out_lw_rad(tmin, tmax, sol_rad, cs_rad, avp)  # rnl
    net_rad = pyeto.net_rad(ni_sw_rad, no_lw_rad)  # rn = rns -rnl
    z = 10
    ws = pyeto.wind_speed_2m(ws_10, z)
    svp = pyeto.svp_from_t(tavg)
    delta_svp = pyeto.delta_svp(tavg)
    # print((net_rad, tavg, ws, svp, avp, delta_svp, psy, shf))
    et0 = pyeto.fao56_penman_monteith(net_rad, tavg, ws, svp, avp, delta_svp, psy, shf)
    # print(et0)
    return et0,ws








rootdir = r'C:\Users\Nical_lmy\PycharmProjects\test\test\8'


rootdir4 = r'C:\Users\Nical_lmy\PycharmProjects\test\test\4'
rootdir5 = r'C:\Users\Nical_lmy\PycharmProjects\test\test\5'
rootdir6 = r'C:\Users\Nical_lmy\PycharmProjects\test\test\6'
rootdir7 = r'C:\Users\Nical_lmy\PycharmProjects\test\test\7'
rootdir8 = r'C:\Users\Nical_lmy\PycharmProjects\test\test\8'



list = os.listdir(rootdir) #列出文件夹下所有的目录与文件

list4 = os.listdir(rootdir4)
list5 = os.listdir(rootdir5)
list6 = os.listdir(rootdir6)
list7 = os.listdir(rootdir7)
list8 = os.listdir(rootdir8)



n = 1
for i in range(0, len(list)):
    path4 = os.path.join(rootdir4, list4[i])
    path5 = os.path.join(rootdir5, list5[i])
    path6 = os.path.join(rootdir6, list6[i])
    path7 = os.path.join(rootdir7, list7[i])
    path8 = os.path.join(rootdir8, list8[i])
    if os.path.isfile(path4 and path5 and path7 and path8):
        # print(path)
        f4 = open(path4, 'r')
        f5 = open(path5, 'r')
        f6 = open(path6, 'r')
        f7 = open(path7, 'r')
        f8 = open(path8, 'r')
        alllines4 = f4.readlines()
        alllines5 = f5.readlines()
        alllines6 = f6.readlines()
        alllines7 = f7.readlines()
        alllines8 = f8.readlines()
        fp = 'E:\\papers_and_notes\\ET_algorithm_and_data\\Shanxi_32766.txt'
        f = open(fp, 'a+')
        # f.write("id s_code year month day tavg tmax tmin ws_10 sunshine_hours rh_mean et0\n")
        for j in range(0, len(alllines4)):
            alllines4[j] = alllines4[j].strip('\n')
            alllines4[j] = alllines4[j].split(' ')
            alllines5[j] = alllines5[j].strip('\n')
            alllines5[j] = alllines5[j].split(' ')
            alllines6[j] = alllines6[j].strip('\n')
            alllines6[j] = alllines6[j].split(' ')
            alllines7[j] = alllines7[j].strip('\n')
            alllines7[j] = alllines7[j].split(' ')
            alllines8[j] = alllines8[j].strip('\n')
            alllines8[j] = alllines8[j].split(' ')
            if alllines7[j][9] == '32766' or alllines7[j][8] == '32766' or alllines7[j][7] == '32766' or alllines8[j][7] == '32766' or alllines5[j][1] == '32766' or alllines5[j][3] == '32766' or alllines5[j][4] == '32766' or alllines5[j][5] == '32766' or alllines5[j][6] == '32766' or alllines5[j][7] == '32766' or alllines6[j][7] == '32766':
                n = n + 1
                continue
            elif alllines4[j][0] =='57411': # 北京站点号
                tmin = float(alllines7[j][9])/10

                tmax = float(alllines7[j][8])/10
                ws_10 = float(alllines8[j][7])/10
                tavg = float(alllines7[j][7])/10
                latitude = math.modf(float(alllines5[j][1])/100)[1]+math.modf(float(alllines5[j][1])/100)[0]/0.6
                altitude = float(alllines5[j][3])/10
                year = int(alllines5[j][4])
                month = int(alllines5[j][5])
                day = int(alllines5[j][6])
                rh_mean = float(alllines5[j][7])
                sunshine_hours = float(alllines6[j][7])/10
                et0,ws = pm(tmin, tmax, ws_10, tavg, latitude, altitude, year, month, day, rh_mean, sunshine_hours)
                f.write(alllines4[j][0]+' '+alllines4[j][4]+' '+alllines4[j][5]+' '+alllines4[j][6]+' '+
                        str(tavg)+' '+str(tmax)+' '+str(tmin)+' '+str(ws)+' '+str(sunshine_hours)+' '+str(rh_mean)+' '+str(round(et0, 3))+'\n')

                print(i, n, j, year, month, day, latitude, altitude, ws_10, tavg, tmax, tmin, sunshine_hours, rh_mean, et0)
                n = n + 1
        f4.close()
        f5.close()
        f7.close()
        f8.close()
        f.close()
